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ABSTRACT 


Hypersonic flow about a body of revolution with or without an 
angle of attack Is of much Interest both from the design and computa- 
tional view points. In this paper a formulation of the complete 
Navier-Stokes problem for a viscous hypersonic flow In general 
cuivilinear coordinates is presented. This formulation is applicable 
to both the axially symmetric and three-dimensional flows past bodies 
of revolution. The equations for the case of zero angle of attack 
have been solved past a circular cylinder with hemispherical caps by 
point SOR finite difference approximation. The free stream Mach 
number and the Reynolds number for the test case are respectively 
22.04 and 168888. The whole algorithm is presented in detail along 
with the preliminary results tor pressure, temperature, density and 
velocity distributions along the stagnation line. 
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1. INTRODUCTION 


The problem of predicting the flowfield about arbitrary blunt 
bodies traveling at supersonic and hypersonic speeds has been the 
subject of much theoretical and experimental research. Supersonic 
flow about blunt bodies is characterized by a detached bow shock wave 
around the nose of the body that divides the flowfield into regions of 
supersonic and subsonic flow. The existence of these regions of mixed 
flow has limited existing analytical methods for the blunt body problem 
to semiemnlr' \il techniques or inviscid flow approximations. Hovever, 
these methods cannot simultaneously predict flow properties in all 
regions of the flowfield. In addition, the relationship between the 
shape of the body and the shape and location of the bow shock has 
restricted most existing solutions to relatively simple body shapes. 

In an effort to overcome the problems associated with existing 
approximate methods, several researchers have developed numerical 
solutions of the full set of governing equations for fully compressible, 
viscous flow, viz., the Navier-Stokes equations, by finite difference tech- 
niques. Unfortunately, the problems inherent in numerical solution of 
the Navier-Stokes equations have usually outweighed the techniques' 
usefulness as a research tool. In particular, the numerical insta- 
bilities that arise from replacing the actual equations with finite 
difference approximations and the large computational fields required for 
accurate resolution of the flowfield have, in the past, prevented the 
treatment of problems of practical interest, i.e., arbitrary body shapes 
and high Reynolds Number flows. However, in recent years, the continuing 
improvement of digital computers and the introduction of more efficient 
numerical methods have made treatment of realistic problems feasible. 


The present Investigation presents a method for the numerical solution 
of the Navier-Stokes equations for both supersonic and hypersonic viscous 
flow about arbitrarily shaped blunt bodies. 

The techniques used in this research are an extension of the 
methods introduced by Warsi, Devarayalu, and Thompson [1,2] and 
Devarayalu [3] for supersonic flow about 2D blunt bodies. The Navier- 
Stokes equations are first transformed from Cartesian coordinates to a 
set of general coordinates by applying basic relations from tensor 
analysis. The numerical techniques developed by Thompson, Thames, and 
Mastin [4,5.6] are then used to generate a boundary fitted curvilinear 
coordinate system for prescribed body and outer boundary shapes. This 
type of coordinate system will map into a rectangular coordinate grid 
in the transformed plane. Since the body and the outer boundary become 
straight lines in the transformed plane, the statement of boundary 
conditions for arbitrary body shapes is greatly simplified. In addition, 
this method of generating coordinates allows coordinate lines to be 
concentrated in any region in the field. A dense mesh system can then 
be imposed in regions of large gradients in the flow variables such 
as boundary layers or across shock waves. 

In the current investigation, a new approach to the transformation 
of the Navier-Stokes equations has been employed. Warsi [7] has shown 
that the transformed equations can be written in a form that retains 
the conservation law or "divergence" form of the original equations for 
all cases except axisymmetric flow. Conservation equations are favored 
in compressible flow problems because they will satisfy the Rankine- 
Hugoniot relations for shock waves when the differencing scheme applied 
to the equations is conservative. 

The next step in the solution process is the development of an 
algorithm to solve the transformed equations numerically. In this 
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research, the solution of the transformed equations was accomplished by 
replacing the partial derivatives in the equations with finite difference 
approximations that produce a set of implicit finite difference equations. 

These equations are then solved by the Successive Over Relaxation (SOR) 

Iterative technique. The implicit formulation was chosen for its uncon- 
ditional stability. The conservation law form of the equations allows 
construction of simple, easily programmed S.O.R. algorithms. 

In the present investigation, the Navier-Stokes equations were 
solved at all points in the field. The bow shock appears in the field 
as a zone of transition of the flow variable smeared over a few grid 
lines. This technique is called "shock capturing." However, several 
authors [8,9,10] have found this technique to be impractical for large 
Reynolds Number flows because of the numerical difficulties associated 
with the large gradients that appear in the flow variables in the regions 
where the bow shock forms. Instead, they treat the shock as a boundary 
across which the Rankine-Hugoniot equations are applied. The Navier- 
Stokes equations are solved only in tre region behind the bow shock. 

The flow in front of thr shock is considered inviscid. This technique 
is called "shock fitting." While generally favored, implementation of 
the "shock fitting" technique can be quite complicated for arbitrary 
body shapes. Therefore, in spite of the numerical difficulties associated 
with shock-capturing, this technique was used in the present research 
because of its simplicity. 

Numerical solutions of the Navier-Stokes equations by finite 
difference techniques are plagued by various types of numerical insta- 
bilities. These instabilities appear as oscillations in the flow 
variables that, if left unchecked, will cause the solution to diverge. 

These oscillations can be controlled by introducing dissipation into 
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the solution by either dissipative differencing schemes or explicit 
artificial viscosity terms added to the calculation of the viscosity 
coefficient. Both of these methods were used in the present research. 

All of the techniques described in this section were used to 
develop a computer program [11] to solve the Navler-Stokes equations 
for the case of axisymmetric flow of a perfect 'ms about a blunt body 
traveling at hypersonic speed. The axisymmetric flow case was chosen 
in order to provide Initial data for an eventual extension of the program 
to 3D flow. The body used in this research consists of a circular 
cylinder with two hemispherical caps. The freestream Mach Number is 
22.04. The Reynolds Number ,>ased on body diameter is 168888, and the 
Knudsen number for this flow is about .0002. 
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2. FORMULATION OF THE PROBLEM 


For the purpose of solving the blunt body problem In general body 
oriented curvilinear coordinate systems, we consider the non-dimensional 
Navier-Stok.es equations in the invariant tensor form. The formulation 
presented here is for a real gas which is assumed to be thermally though 
not necessarily calorlcally perfect. A computer program based on these 
general considerations has been developed which can also be easily used 
tor calorlcally perfect gases. 

A dimensional quantity Is denoted by a superscript * and the value 
of a quantity at free stream by a subscript °°. The equations have been 
non-dlmensionallzed by using the following free stream values. 

* 

L for all lengths 

* 

for velocity vector 

A 

density 

* * 2 

p U for pressure and energy per unit volume 

00 00 

* 

T^ temperature 

A 

viscosity 

* 

C specific heat at constant pressure 

The non-dimensional form of the conservation equations are 

( 2 . 1 ) 
(•■*.. 2 ) 
(2.3) 


+ div(pv) - 0 
|^(pv) + div i-O 
U + div b - 0 


J 


— 


. . 


where 


t ■ pvv + pi - to 


c ■ KI + d 


K - X div y 


d ■ p def v 


b ■ (f + p)v - eo*v - Sp grad T 


* - pe + y p|v| 


* , * * * 

e - 1/R e - ^/P.U.L 


* * . * 2 

S - eTCC/PU 

00 poo p r OD 


* * * 

P ■ |i C /k 
r P 


(2. A) 


k - dimensional fluid conductivity. 


In Eqs. (2.4), p and X are the non-dimensional first and second coefficients 
of viscosity respectively. 

The set of equations (2.1) - (2.4) form a closed system when the 
following additional conditions are associated with them. 

(i) The equation of state 

p - C 1 (y - l)pTC p /y (2.5) 


C 


1 


* * 

C T /U 

nao OD 


*2 


where 


Cp ■ non-dimensional specific heat at constant preasure 


y - ratio of specific heats . 


Both Cp and y ere assumed to be functions of temperature T. 


(ii) Sutherland's viscosity law: 


p - (1 + S 1 )T 3 / 2 /(I + S^ 


( 2 . 6 ) 


where S. - D*/T* , D* * 110.33°K . 
1 00 


(ill) Stokes' law: 


3A + 2 p ■ 0 


(2.7) 


(iv) Eucken's formula: 


P r (T) - A y /(9 y - 5) 


( 2 . 8 ) 


The temperature T is obtained by solving a fifth degree equation defined as, 
(see appendix). 

(A x - 1)T + y(A 2 T*)T 2 + y(A 3 T* 2 )T 3 


1 *3 . i *4 

+ A< A 4 T oo >* + j(A 5 T )T* 




(2.9) 


where 


c 2 - c 1 /f(i) 


f (T) - A x + (A 2 TJT + (A 3 T*‘)T 2 


+ (A 4 T* 3 )T 3 + (A 5 T*V 


( 2 . 10 ) 
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II IU II * an ■ IJ | 




— 


The conatante are 

A 1 - 3.6535 r ' 

Aj - -1.33736 x 10" 3 /°K 
A 3 - 3.29421 x 10~ 6 /(°K) 2 
A 4 - -1.91142 x 10" 9 /(OK) 3 
A 5 - 0.275462 x 10“ 12 /(°K) 4 

On transformation, Eqs. (2.1) - (2.3) have been expressed in the 
strong conservation-lav form (cf. Uarsl [7]) as follows 


3o . 3 , 1. 

3? + -T (ov > 
3 C 


( 2 . 11 ) 


^ < ° Y> + 7^ (xik?1> " 0 


f + ii;. o 

3t 3C 1 


(2.12) 


(2.13) 


Where £* (1 ■ 1,2,3) are the curvilinear coordinates, v* are the contra- 
variant components of the velocity vector v, and a^ are the covaviant base 
vectors. Further, 


o - /g p 

X ik - ovSr* + g ik (P - ex) - cD ik 

Y 1 - (E + P)v l - c* 1 - Sv<g g 1J 2L- 

K 2 
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(2.14) 


P-« / gP.E"* / g'l , .X*/glC 

• ^glKv 1 + w(g ik g. v v 1 )v" 
ja ,k ,m 

D 1 * - Si d 1 * 


and g la Che determinant of the coefficient matrix g^. 

2.1 Axlally-Symstetric Body 

We now particularise the equations for the case of flow past an 
axially syometric body of arbitrary shape. Referring to Fig. 1 let 
■ C, £ 2 ■ n be the coordinates in the meridian plane and £ 3 E <> the 
azimuthal angle, where x is the axis of symmetry. Thus the generation 
of the coordinates ( and n, with n ■ constant along the body contour, is 
exactly a two-dimensional problem. The forms of and (the 

Christoff el symbols) for the axially symmetric body ate: 

o ■ x 2 r 2 

g ll C C 


l> m 

g 22 n n 


F 33 - r 2 - y 2 + t 2 


J ■ x.r - x r. 
C n n C 


g 12 " Vn * Vn 

g 13 " g 23 " ° 

g - r 2 lBn 822 " " r2j2 

g11 * r2 *22 /g 
g2 ' r2g ll /g 
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g 33 . i/ r l 


g u . -r 2 g 12 /g 


.13 m „2 3 B 0 


t r 2 3> n ,, !!ii _ !!u 

11 " 2g g 22 3C " g 12 ( ‘ 3£ 3n 


)] 


12 " Ti la 22 an* ' 8 12 H 




r 1 

* 1 1 


li 


3 g 


22 , 


r 1 P 1 ■ 0 

*13 23 


r 1 

*22 


_2 3 g. 2 ^*22 

2g lg 22 (2 ~3n aT _) 


3g 22 1 
*12 3n 


r 33 " r (g 12 r n “ g 22 r C ) 


r il " 2g lg U (2 3C 


Sg 12 3g ll 


3g 


11 , 


3n * “ g l2 3C 1 


.2 *g r ' g 12 Uft 22 , 


12 * £ lg ll ^ - *12 (2 -->T )] 


r 2 3 * 22 

r i2 “ 2g (g ll'H ' 8 12 3n 


13 ‘23 


r 2 - 0 


r 33 " r (g !2 r 4 • g H r n > 


r 3 • r 3 ■ r 3 ■ tL ■ o 
11 *22 12 33 


r 


3 

13 



r 
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r 3 m — n 
*23 r 


(2.15) 


The base vectors are 


?! * **£ + Jr^cos** + kr^sinij) 


&2 * + Jr^coscfi + kr^sinifr 


(2.16) 


a^ ■ r(kcos$ - jsin$) 


where i, j, k are the unit vectors along the Cartesian coordinates x, y, 
z, respectively. 

The relations between the Cartesian components of the velocity vector 
denoted as 0, V, U and the contravariant components are 


U = v*x, + v 2 x 


V * Scos<|> - rv 3 sin<t> 


W ■ Ssin$ + rv 3 cos$ 


(2.17) 


where 


S ■ v*r„ + v 2 r 


Vcos<|> + Wsin$ 


From (2.17), we have 


(2.18) 


v 1 = (Ur - Sx )/« / g/r 2 

n n 
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(2.19) 


v 2 - (Sx^ - Ur^)//g/r 2 
v 3 ■ (Wcos«ji - Vsin<|i)/r 

It must be noted that rv 3 is the velocity of swirl. 

ik 

The symmetric momentum- stress tensor X appearing in F.q. (2.12) 

ik 

and defined in (2.14) is now used to define another tensor XM defined 
below. 

XM 11 - X n x r + X 12 x 

4 n 

XM 12 - X 12 x + X 22 x 

4 n 

23 

XM 13 - X 13 x r + X x 

4 n 

XM 21 = (X ll r f + X 1,: r t )cos(|i - rX l3 sin<{> 

XM 22 = (X 12 r f + X 22 r (i )cosi|i - rX 23 sini}> 

XM 23 * (X 13 r f + X 23 i^)cos,}> - rX 33 sin<j> 

XM 31 * (X n r f +■ X^r^Jsini)) + rX* 3 cos<J> 

XM 32 = (X 12 r. + X 2- r )sin4> + rlv 3 cos$ 

4 n 

XM 33 = (X 13 r. + X 23 r^)sini> + rX 33 cos4> (2.20) 

When Eqs. (2.1b), (2.17) and (2.20) are used in Eqs. (2.11) - (2.13), 
we get the following equations. 

I? + H (t,vl) + ^ (ov,!) + li (ov ' ) ‘ 3 <2 - 21) 
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, O 3„ / dv 3 3v^ 

+ f u - £i,(x t ir-'.ir 1 


. sHi {r isi + r 5ii) . o 

r 2 vr g 3$ n 3<j> ' 


(2 . 2 / ») 


3E . 3Y 1 . 3Y 2 3 r , ,,3v 3 . i 

at IF '5iT “ E F 1,rJ( aT + ? )v 


t 3 3v1 ur 3 v 3 , 3v J 3v J ., 

+ UrJv 34> + J 8 22 3£ 8 12 3n 


3 r-t w3v 3 , ? . „ 39v 2 xl 

E F l>rJ( aT + 7 >v + 1,rjv V >’ 


+ H. 3 ''2 (g i*i _ . iyi) i 

J ^ g ll 3n 8 12 35 >] 


+ ^[r(E + P)v 3 -e* 3 -^ff]=° 


(2.25) 


In numerical vector form, Equations (2.21) - (2.25) can be written 


as 


3W 

at 


+ 


K 

8C 



3<f> 


0 


(2.26) 


For the case of axially symmetric flow at zero angle of attack and 
with zero swirl the preceding equations become simplified by setting v 3 = 0, 
and ft = 0 


Boundary Conditions : 

(i) Free stream conditions: 

The flow at upstream infinity is assumed to be uniform at supersonic 

* -k * 

free stieam conditions U , p , p at an angle of attack a with the axial 

ao oo a? 

direction. The non-dimensional form of the free stream conditions are 
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p 


1 


p * p 
1 00 


M - M 

QD 


U ■ cosa , S ■ sinacosif , rv 1 ■ -sinasin<f> 

(ii) At the body surface: 


(2.27) 


U - S - v 3 = 0 


T 


T or 
w 


(i>I) 


(2.28) 


(iii) Out-flow plane: zero derivative conditions. 

(iv) Symmetry conditions on the planes $ » 0, $ ■ u for the axially 
symmetric case. 

For the case of aero angle of attack, the calculations on the stagna- 
tion line £ = £ are to be performed separately. Let F be either U, p or 
s 

E. Then on r he stagnation line 


FU s ,n,$) - FU s ,n,<M-it) 


(II) „ _ ( il) 
ac t-*i 8 




(J> = lf+Tl 


(2.29) 
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3. NUMERICAL SOLUTION OF THE NAVIER-STOKES EQUATIONS 


The numerical methods used to solve the Navier-Stokes equations 
are presented in this section. The development of finite difference 
approximations of Equations (2.26) and the algorithm used to solve 
them are described. Boundary conditions and the calculation of values 
along the stagnation lines are discussed. The procedures used to 
control nonlinear instabilities in the solution are described. 

Finally, a discussion of the initial conditions for the test case used 
in this research is given. 


A. Finite Difference Approximations and the Solution Algorithm 

The numerical vector Equation (2.26) for the case of zero angle of 
9E 

attack, — = 0, was discretized into a set of difference equations by 
replacing the partial derivatives with finite difference approximations. 
These difference equations were written in a form that allows either of 
two fully implicit differencing schemes to be used. In both of these 
schemes, spatial derivatives were replaced by second order central 
difference approximations. In the first scheme, the time derivative is 
replaced by a three point backwards difference approximation at N+l 
where N is the time step of size At. 

In the second scheme, the time derivative is replaced by a first 
order backwards difference at N+l. The spatial derivatives are 
replaced by the average of the derivatives evaluated at time steps 
N+l and N. This scheme is commonly referred to as "trapezoidal" or 
Crank-Nicolson differencing. Both of these differencing schemes are 
second order accurate m space and time. 

SUR iteration was used to solve the system of equations resulting 
from the difference approximations described above. The SOR algorithm 
for both differencing schemes is 


lo 


(3.1) 


£i - C p) ♦ “V 


where u> is the relaxation parameter and (p) denotes values at the 
previous iteration. The residual vector, R, is 




(i+6 3<j - - «?! 1 , (p) - ‘2“ 81,5 


N+l 

I.J 


" (6 2- 6 3 )RHS I,J 


(3.2) 


where 


^I.J " -I+l.J ' -I-1,J + -I.J+l " >I,J-1 _2 -I,J (3.3) 

In equation (3.2), setting the parameters, ^2 = *^3 = 3/3, yields the 
three point backwards time differencing scheme. Setting, 6^ = 0 and 
&2 = 1/A, yields the trapezoidal scheme. 

All parameters in Equation (3.2) evaluated at N+l use the most 
recent values of W from the previous iteration. Modifications to 
Equation 0.2) to control non-lineav instabilities are discussed later 
in this chapter. 

B . Boundary Conditions on t h e Infinity Boundary 

The infinity boundary is divided into inflow and outflow region?. 

The freestream values of the flow variables are maintained along the 
inflow boundary. Values along the outflow boundary are computed by 
extrapolation of the physical flow parameters (p, U, S, V) from the 
field. The extrapolation formulae for the outflow boundary is derived by 
assuming the spatial derivatives of the flow variables along lines of 
constant £ at the outer boundary are equal to zero. The three point 
backwards difference approximation of the first derivative then can 
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be used to write the following formula: 


F I,JMAX " * 4F I,JMAX-1 " F I,JMAX-2 )/3 (3,4 * 

where F ■ (p, U, S, Y) and J ■ JMAX denotes the outer boundary. 

C. Boundary Conditions Along the Body Surface 

The body surface is taken to be an impermeable, isothermal boundary. 
The no slip condition (U » S ■ v 1 ■ v 2 ■ 0) is imposed along the boundary. 
Pressure and energy at the boundary are computed by Eqs. (2.5) and (2.9). 
The wall density is computed by solving the continuity equation at the 
body surface. The no slip condition permits the continuity equation at 
the boundary to be written as 


(l£.) . _ (iiL) 

'8c 1,1 v 9rr 1,1 


(3.5) 


where B ■ ov 2 and J ■ 1 denotes the body surface. 

Two different differencing techniques were used to evaluate 
Equation (3.5). In the first technique, the time derivative is evaluated 
at the N+l time level by a first order backwards difference. The spatial 
derivative is evaluated at the Nth time level by a three point backwards 
difference approx imat ion. Equation (3.5) then becomes 


N+l N At .,_N „N , 
’l.l ■ °I,1 ■ r l4B I,2 ' B I,3> 


(3.6) 


N 

where b” ■ 0. This equation was used for the first two time steps of 
the solution. Equation (3.6) is an explicit equation and is known to 
cause instabilities with repeated use. Therefore, an implicit equation, 
derived by replacing both the time and the space derivatives with three 
n^nt backwards differences evaluated at time level N+l, was used at 
succeeding time steps. These approximations yield 
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f 

I 

; 

i 




i 


1 


N+l 1 N N-lv 
°I,1 M 3 (4 °I,1 " °I,1 


At ( , N+l _ N+l. 
3 (AB I,2 B I ,3 


(3.7) 


where B 


N+l 

1,1 


D. Calculation of Values Along the Stagnation Lines 

For axisymmetric flow, Equations (2.26) contain singularities along 
the forward and rear stagnation lines where r • 0. The values of the 
flow variables along the stagnation lines must therefore be computed 
by extrapolation from the field. The extrapolation scheme proposed by 
Windhopf and Victoria [12] was used in the present research. Let F(£,n,$) 
be any one of the flow parameters p, U, or Y on the stagnation lines. 

Then, by symmetry 


F(I .J,*) - F(I ,J,$+tt) 

b S 


(3.8) 


and 


(^) 


- <^) 

3C <f>*4>+ir 


(3.9) 


where I g a 1 at the stagnation line. 

An extrapolation formula for the flow variables along the stagnation 
line can be constructed by replacing the derivatives in Eq. (3.9) by 
three point forward difference approximations. This substitution and 
Eq. (3.8) combine to yield 

F I ’ft I4(F T + F I ) ~ (F t + F, )] (3.10) 

s.J s+l,J s-l,J s+2,J s-2,J 


If the symmetry of the flow variables is ideally maintained during the 
course of the solution, Eq. (3.10) can be replaced by 


V J 


5 WF W 


F I j) 

s+2 J 


(3.11) 


1 
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Equation (3.11) was initially used to compute the values of the flow 
variables along the forward stagnation line. 


Along the rear stagnation line, the values of the flow variables 
were computed by averging the values at I"2 and IMAX-1. The Cartesian 
velocity component, S, is set to zero along both stagnation lines. The 
contravariant components of velocity, v 1 and v 2 , car. then be computed 
by Equations (2.19). For S ■ 0, Eqs. (2.19) yield 


v 


1 



v 2 



J 


(3.12) 


E. Control of Numerical Instabilities 

The numerical instabilities described in Chapter 1 are an inherent 
part of numerical solutions of systems of partial differential equations 
such as the Navier-Stokes equations. The oscillations in the flow 
variables that signal the onset of numerical instabilities must be 
damped as soon as they appear in the solution in order to maintain a 
stable solution. Although several methods for controlling numerical 
instabilities have been proposed, the most successful techniques for 
compressible flow problems have introduced terras into the finite dif- 
ference approximations of the equations being solved to "filter" or 
"smooth" the oscillations. These additional terms act to add artificial 
viscosity that increase the effects of dissipation in the solution. An 
equivalent method is to add an explicit artificial viscosity term to the 
calculation of the coefficient of viscosity, u. Both of these methods 
were used in the present research. 

Initially, the Shuman filtering technique employed by Vi iegenthart 
[13] and Marten and Zwas [14] was used. The Shuman filter is implemented 
by replacing and * in Eq. (J.2) with 




J 
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w 

"i.j 


+ 2. r.J 


*-I+l,J + -I-1,J + -I t J+l + 


(3.13) 


where j represents the physical values of the flow variables 
(p,pU,pS,¥) and a is a switching parameter that can be used to turn 
the filter on and off or to vary the strength of dissipation. Setting 
a equal to 0.5 yields the form of the Shuman filter used in Kefs. 

[1,2,3]. Application of this form of the filter at all points in the 
field will eliminate oscillations in the flow variables. However, 
the excessive dissipation Introduced by the full Shuman filter lowers 
the effective Reynolds number of the flow and delays convergence to a 
steady state. Therefore, it is desirable to apply the filter only in 
regions where oscillations occur. This can be accomplished by setting 
the value of the switching parameter by either of two methods. 

In the first method, the field is searched prior to the start of 
a time step for N shaped wave forms. The filter is then applied at 
the inner points of the N-wave. The second method is similar to the 
technique used in Ref. [ 14 ] . The value of a is varied from point to 
point in proportion to the magnitude of the local gradient of one of 
the flow variables such as density or pressure. Therefore, more 
diffusion will be applied in regions of high gradients, such as across 
shock waves, where numerical instabilities have a tendency to occur. 

It should be pointed out that the N-wave and the switched forms of the 
Shuman filter will not completely eliminate oscillations enough to allow 
the solution to converge the steady state. 

Because of the severity of the full Shuman filter, a second technique 
for introducing dissipation into the solution was tried. This technique 
adds an explicit fourth order artificial viscosity term to the computation 
of the viscosity coefficient. This method of adding artificial viscosity 
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is related to the fourth order smoother used by Baldwin and MacCormack 
[15] and was used successfully by Thompson [16] in ar incompressible 
flow problem. In the current research, this method was implemented by 
replacing the viscosity computed by Sutherland's law with 

V - u U ♦ C(— ♦ — )] <3.1ij 

H 2 an 2 

where P s is the viscosity given by Eq. (2.6). C is an arbitrary constant 
whose magnitude is determined experimentally. 

F. Initial Conditions 

The uniform flow conditions used in this research were as follows: 

M ■ 22.04, |v I ■ 7000.0 m/sec, p ■ 7.73067 * 10 kg/m 3 , and 

cd * u 1 * On 

A 

T^ * 250.8l°K. From these conditions the pressure, coefficient of 

viscosity and Knudsen number are computed to be: P v 5.57043 N/m* , 

p * 1.60209 « 1C kg/m-ser and K ■ 0.000193. For a body diameter 
“ n 

of five meters, the free stream Reynolds number is equal to 168888. 

The ratio of specific heats was assumed to be y ■ 1.41. An isothermal 

* 

wall temperature, T y ■ 1000. 0°K, was maintained throughout the course 
of this research. 

Because the solution algorithm used in the present investigation 
is iterative, an initial guess for the values of the flow variables 
is needed at all field points in order to start the solution. However, 
finding a method of starting the solution that didn't diverge in the 
initial time step proved to be one of the major tasks of this research 
project. Although several methods of starting the solution were 
attempted, only one technique provided an initial guess that would 
converge for the first time step. In this technique, the U component 
of velocity is assumed to vary linearly along a line of constant f. 
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from the body to the outer boundary. The S component of velocity was 
assumed to be zero. Temperature, pressure and density at all points In 
the field was assumed to be equal to their values along the outer 
boundary. A similar technique was applied in Refs. [8,9,10] to ini- 
tialize the values of the flow variables in the region between the 
body and the bow shock wave. In the present research, it was found 
that this method would converge if the time step size was lowered to 
At ■ .001. The acceleration parameter was 0.9 for all four equations 
of motion. Although the first time step converged using this method 
of starting, the fulJ Shuman filter was required for about ten time 
steps to maintain stability 
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k. COMPUT AT I ON AI. PROCEDURES AND RESULTS 
k. Computational Procedure s 

The numerical techniques described in the previous chapter were 
U3ed to write a computer program that can solve the Navi er-Stok.es 
equations for either 2D or axisymmetric flow about an arbitrary axially- 
synunetric blunt bodv. The solution process begins with the generation 
of a coordinate system and the calculation of the corresponding metric 
data for a prescribed body and an outer boundary cf. Fig. 2. As shown 
in Fig. 3, the outer boundary used in the present research is an ellipse 
centered at the midpoint of the body. It is assumed that the outer 
boundary is placed far enough away from the body surface so that the 
flow on the downstream boundary will be supersonic. The region between 
the sixth and the forty-sixth lines of constant £ is taken as the outflow 
boundary. 

Once the metric data for a given coordinate system is computed 
and stored on file, the actual solution of the Navier-Stokes equations 
is begun by assuming an initial guess of the solution for the entire 
computational domain. The solution is then advanced through time until 
a steady state solution is reached. However, as has been previously 
mentioned, finding an initial guess that did not diverge in the first 
time step proved to be one of the major tasks of this research. 

It was originally intended to first obtain a steady state solution 
for the case of 2D flow and then use the resulting data as an initial 
guess for the axisymmetric flow case. Cenerally, the uniform flow 
conditions along the upstream boundary an be used to provide an 
initial value at each field point since the initial guess and the 
resulting transient solution are not required to be physically 
realistic. For this initial gucs.s, the boundary conditions along the 

2U 




wall become the driving functions for the solution. However, as in 
Ref. [3], it was found that this method would not work when the free 
stream Mach number is very high. After considerable numerical experi- 
mentation, it was found that the solution could be started at the 
desired Mach number (22.04) by assuming a linear variation in velocity 
between the body and outer boundary along lines of constant £. Energy 
and density were set to the free stream values at all points in the 
field. However, this technique had two disadvantages. The first dis- 
advantage was that the time step size had to be lowered to .001 for a 
fixed acceleration parameter of 0.9 for all the four equations of 
motion. Extensive numerical experimentation did not produce a combi- 
nation of time step and acceleration parameter that would allow a 
larger time step to be used. The second disadvantage was that the 
values of temperature, velocity, and pressure in the region between 
the expected location of the bow shock and the outer boundary produced 
by the linear velocity distribution differed substantially from the 
free stream values. The formation of a shock transition region was 
therefore substantially delayed because of the time required for the 
free stream flow to convect into this region and increase the magnitude 
of the velocity toward the free stream value at a time step of .001. 

The 2D case was started using the linear velocity distribution and run 
for several hundred time steps. Although the solution showed signs of 
eventually approaching the free stream in the region in front of the 
shock, the computer time required to obtain a steady state solution 
would have been considerable. 

Figures 4 and 5 present distributions of temperature and velocity 
magnitude along the stagnation line that are typical of the linear 
velocity distribution initial guess. 
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At this point in time, the cost of using the local computing 
facility became prohibitive. Further progress was delayed until 
arrangements could be made to use the computational facilities at 
NASA's Marshall Space Flight Center. Because of the delay induced by 
the transfer of the computer program to NASA-Marshall and the time 
limits Imposed by contractual requirements, it was decided to abandon 
attempts to obtain a steady state 2D solution. Instead, the methods 
described above were used to restart the solution for the case of 
axisymmetric flow. In an attempt to speed up the solution process, 
a transition region between the values of the flow variables at the 
anticipated location of the bow shock and the free stream values was 
introduced explicitly after fifty time steps. The stand off distance 
of the shock along the stagnation line was computed using the empirical 
equations of Billig [17]. A linear variation between the values of 
the flow variables at the n location directly in front of the anticipated 
shock location and the free stream values was taken over five cell 
lengths in the q direction. The solution was continued in a normal 
manner from this point on. 

It was found that the full Shuman filter described in the previous 
chapter was required to control instabilities induced in solving the 
finite difference equations. However, once the flow field was established, 
it was found that repetitive use of the full Shuman filter tended to 
smear the shock transition region toward the outer boundary. Unfortunately, 
t he dissipation introduced by the full filter was required to control the 
"wiggles" that appeared in the flow variables. The switched and N-wave 
forms of the filter did not appear to introduce enough dissipation to 
control these oscillations. 

At this point, the calculation of the explicit artificial viscosity 






term described in the previous chapter was introduced into the program. 

It was found that the arbitrary constant, C, had to be on the order 
-6 

of 10 for the method to be effective. This method of introducing 
dissipation seemed to work quite well at first. However, as the 
solution progressed in time, the value of the temperature along the 
stagnation line in the region where the shock wave had started to form 
suddenly started to decrease. This trend continued until a negative 
value of temperature was reached. It is felt that modifying the 
viscosity coefficient in this manner produced unrealistic values of 
energy that lead to the calculation of the negative temperature. 

The solution has been run through 700 time steps at a step size 
of 0.001. Efforts to increase this step size have been unsuccessful. 
For a fixed acceleration parameter of 0.9, it took on the average 
.825 min. of computer time to converge a time step to a tolerance of 
10 J . Convergence was usually obtained after two or three iterations. 
For a 51 x 50 mesh, this time represents a computational effort of 
about .0194 sec. per grid point per time step. 

All runs were made on UNIVAC 1100/80 series computers at 
Mississippi State University and NASA-Marshall . The computer program 
was written in FORT? AN using the UNIVAC ASCII FORTRAN compiler. This 
program required approximately 102,000 words of core on the 1100/80. 

b. Discussion of Results 

The solution of the Navier-Stokts equations for the case of 
axisymmetric flow with the initial conditions given in the previous 
chapter has been obtained through a non-dimensional time of 0.7. 

Since this is still a very early time, the results at this time can 
only be considered as preliminary. Therefore, only qualitative 
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judgments can be made about the following results. Figures 6-11 
depicts the non-dimension il values of density pressure, temperature, 
and velocity magnitude along the forward stagnation line (I * 26) 
at a non-dimensional time of 0.7. The trends in these values seem 
satisfactory. The transition region that denotes the existence of 
the shock wave is readily seen. The standoff distance of the shock 
as given by the plots of temperature and pressure compare favorably 
with the standoff distance computed by the empirical equations given 
in Ref. [17]. Density and velocity magnitude are plotted as functions 
of both che non-dimensional distance from the nose of the body and 
the coordinate n- The pronounced rise in density in the viscous 
region close to the body is clearly seen in Figures 6 and 7. 

Figures 8 and 9 show that the velocity profile in the transition 
region has been smeared over about twelve grid cells. The temperature 
profile along the stagnation line is given in Figure 10. The maximum 
value of temperature shown in the plot seems unrealistically high. 

This is thought to be due to a combination of the transient nature of 
the flow and the assumption of an ideal gas. Figure 11 presents the 
distribution of pressure along the stagnation line. Figur^ 12 represents 
a plot of density contours about the body whose values are greater 
than that of free stream values. This plot indicates a large region 
of rarefied flow extending over much of the field. This rarefaction 
is felt to be a transient phenomena induced by the initial guess. 
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5. CONCLUSIONS AND RECOMENDAT IONS 


A method for solving the Navier-Stokes equation for hypersonic 
flow about axially symmetric bodies has been presented. Although only 
preliminary results have been obtained, the following conclusions 
can be drawn. 

(1) The technique offers significant computational advantages 
because of the conservation law form of the equations and 
the reduced amount of metric data required for axisymmetric 
and 3D flow calculations. 

(2) The major difficulty encountered with the method was finding 
an initial guess that would remain stable and converge to 

a realistic solution. Although the linear velocity variation 
technique used in the present research supplied a stable 
initial solution, large amounts of dissipation were required 
to maintain stability. The time step imposed by this 
starting method also is very costly in terms of computer 
time. The method of starting from a steady state 2D 
solution needs to be pursued. A gradual start in which the 
solution is started at a lower Mach number and gradually 
brought up to the desired Mach number should be investigated. 

(3) The problem of controlling numerical instabilities was not 
satisfactorily reconciled. Although the switched form of 
the Shuman filter was not particularly effective in this 
research, it seems to offer the most promise for reducing 

the severity of numerical instabilities without substantially 
retarding the progress of the flow toward steady state. 
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(4) A more refined mesh system should be employed to increase 

the accuracy of the solution in the region of the bow shock. 
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Figure 6. Density Distribution Along the 
Stagnation Line at Tirse 0.7. 
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Velocity Distribution Along the 
Stagnation Line at Time 0.7. 
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Appendix 

Thermodynamic Equations 


On the raolal basis, the specific heat at constant pressure of a 
thermally perfect gas in the range 300°K to 1000°K is given by, [18], 

— * , *3 *4 

C = R(A, + A«T + A_T* 2 + A. T + A C T ) , (A-l) 

p Li J 4 j 

where the constant A^^ have been giv»n after Eq. (2.10), and R is the 
universal gas constant, 

R - 8314.3J/kgraol - °K . 


Introducing 


P m 


R = — 


where m is the molecular weight, we have 


* * *2 *3 *4 

C = R(A. + A-T + A-T + A.T + A,T ) 

p 1 2 3 4 5 


On non-dimensionalization 


. c 


we obtain 


c -iili. 

p f (i) 


where f(T) has been defined in Eq. (2.10) 


(A-2) 


(A-3) 
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Note that 


C - Rf (1) 

poo 


(A-4) 


* A 

The ratio of the specific heats and C denoted as 

p v 

* < * 

Y - C /C 
P v 


can now be expressed a function of T by using the thermodynamic equation 


* * 

C - C - R 
P v 


and the Eqs. (A-3) and (A-4) as 


y(t) 


f(T) 


f(T) - 1 


(A-5) 


For a thermally perfect gas 


a a a 
de = C dT 
v 


(A-6) 


and 


* A A 
dh = C dT 
P 


(A-7) 


On non-dimensionalization, Eq. (A-6) becomes 


de 


A A 

T C C 

oo n 00 p 

— dr- ul 


(A-8) 


Using (A-3) , (A-4) and (A-5) , we get 

de = C 2 [ f (T) - l]dT 
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where, 


c 2 - C^fd) 


C 


1 


RT 

OP 

-S7 


U 


(A-9) 


Substituting (2.10) for f(T) in Eq. (A-0) and integrating, we obtain 
e(T) - C 2 1 (A 1 - 1)T +y(A 2 T*)T 2 + i(A 3 T*^)T 3 

+ 7(A,T* 3 )T 4 + -kA,T* t, )T'] (A-10) 

4 4 <*> 5 5 00 

Having expressed e as a function of T, we now use the expression for e 
from (2.4) to have 


Substituting e(T) 
for obtaining the 



from (A-10) in (A-ll) we get a fifth degree equation 
temperature, viz., Equation (2.9). 
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